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Abstract 

We analyze the timing of photons observed by the MAGIC telescope during a flare of the 
active galactic nucleus Mkn 501 for a possible correlation with energy, as suggested by 
some models of quantum gravity (QG), which predict a vacuum refractive index - 1 + 
(E/MqcriT, n = 1,2. Parametrizing the delay between y-rays of different energies as At - 
±TiE or At = ±TqE^, we find t/ = (0.030±0.012) s/GeV at the 2.5-o- level, and = (3.71 ± 
2.57) X 10"^ s/GeV'^, respectively. We use these results to establish lower limits Mqgi > 
0.21 X 10'^ GeV and Mqg2 > 0.26 x lO" GeV at the 95% C.L. Monte Carlo studies confirm 
the MAGIC sensitivity to propagation effects at these levels. Thermal plasma effects in the 
source are negligible, but we cannot exclude the importance of some other source effect. 



Key words: Gamma-Ray Sources (Individual: Mkn 501), Active Galactic Nuclei, 
Quantum Gravity, Supersymmetric Models, Particle Acceleration 
PACS: 98.70.Rz 03.30.+p 04.60.-m 95.85.Pw 



1 Introduction 



It is widely speculated that space-time is a dynamical medium, subject to quantum- 
gravitational (QG) effects that cause space-time to fluctuate on the Planck time 
and distance scales [1-8], for reviews see [9]. It has also been suggested that this 
'foaming' of space-time might be reflected in modifications of the propagation of 
energetic particles, namely dispersive effects due to a non-trivial refractive index 
induced by the QG fluctuations in the space-time foam. There are microscopic 
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string-inspired models [1-3] that predict only subluminal refraction, and only for 
photons [10], suppressed either linearly or quadratically by some QG mass scale: 



Ac E Ac 



, or — = —. (1) 

MoGl c M2 



One might guess that the scale Mqgi or Mqg2 would be related to Mp, where 
Mp = 2.4 X 10'^ GeV is the reduced Planck mass, but smaller values might be 
possible in some string theories [2, 3], or models with large extra dimensions [11]. 
Superluminal modes and birefringence effects are also allowed in some other mod- 
els [4-8]. 

A favored way to search for such a non-trivial dispersion relation is to compare the 
arrival times of photons of different energies arriving on Earth from pulses of dis- 
tant astrophysical sources [1, 12]. The greatest sensitivities may be expected from 
sources with short pulses, at large distances or redshifts z, of photons observed 
over a large range of energies. In the past, studies have been made of emissions 
from pulsars [13], y-ray bursts (GRBs) [1, 11, 12, 14-17] and active galactic nu- 
clei (AGNs) [18, 19]. In particular, a combined analysis of many GRBs at different 
redshifts made possible some separation between energy- and source-dependent ef- 
fects, and yielded a robust lower limit Mqgi > 0.9 x 10'^ GeV [15]. Astrophysical 
sources that produce very high energy photons in the TeV range or higher could 
improve significantly the sensitivity to Lorentz violation, if one could distinguish 
source and propagation effects. Flaring AGNs are celestial objects with the desired 
properties, and a pioneering study of a flare of the AGN Mkn 421 yielded a sensi- 
tivity to Mqgi ~ 4 X 10^^ GeV [18]0 

In this Letter we analyze two flares of Mkn 501 (z = 0.034) observed by the Ma- 
jor Atmospheric Gamma-ray Imaging Cerenkov (MAGIC) telescope [22] between 
May and July 2005. After applying standard quality checks, data covering a total 
observation time of 31.6 h spread over 24 nights survived and were analysed [23]. 
The data were taken at zenith angles of 10° -30°, resulting in an energy threshold 
(defined as the peak of the differential event-rate spectrum after cuts) of « 150 GeV. 
The air-shower events were subjected to the standard MAGIC analysis [24], which 
rejects about 99% of hadronic background events, while retaining 50-60% of the 
y-ray induced showers. The y-ray energies are, in a first approximation, propor- 
tional to the total amount of light recorded in the shower images; corrections are 
applied according to further image parameters [25] obtained from the analysis. The 
achieved energy resolution is ~ 25% over the range 170 GeV to 10 TeV. The arrival 
time of each event is obtained with sub-ms precision. 



Stronger limits hold in models predicting birefringence [20,21], but these do not apply to 
stringy models of QG-induced vacuum dispersion [2, 3], in which birefringence is absent. 
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During the observations, variations in the y-ray flux by an order of magnitude were 
observed, with the maximum integrated flux above « 150 GeV exceeding (1 1.0 + 
0.3)x 10"^° cm~^s~' . In the two nights with the highest flux, high-intensity outbursts 
of short duration (flares) were recorded with characteristic rise and fall times of 
1 -2 minutes. While the flare of July 9 was clearly visible over the full energy range 
0. 1 5 - 1 TeV and reached a peak flux more than a factor two higher than before and 
after the flare, that seen on June 30 was concentrated at low energies (0.25-1 TeV) 
and less significant. In the analysis below, we applied cuts on the image parameter 
ALPHA [25], describing the gamma shower arrival direction: \ALPHA\ < 10°, and 
on energy: Ey > 150 GeV. 



2 Timing analysis 

The spectral time properties of the most intense portions of the flares were quan- 
tified in [23] using four diff'erent energy bands with boundaries at 0.15, 0.25, 0.6 
and 1.2 TeV, the fourth band extending to infinite energies. In the June 30 flare a 
signal above a uniform background appeared only in the energy band of 0.25-0.6 
TeV, which did not permit any conclusion on the time- spectral properties of the 
signal. For the flare of July 9, a time lag of about 4 minutes was found for the 
maximum of the time profile envelope for photons in the 1.2-10 TeV energy band 
relative to those in the range 0.25-0.6 TeV. The diff'erence between the mean en- 
ergies in these two bands is « 2 TeV, which would lead to a naive estimate of a 
time delay of about 0.12 s for each GeV of energy diff'erence. However, this ap- 
proach is too simplistic, since the energy range covered by the 1.2-10 TeV band 
is much larger than the energy diff'erence between the two bands, so the binned 
estimator used in [23] is inadequate for constraining possible linear or quadratic 
energy dependences. In view of this and the limited number of photons, we im- 
prove here on the binned estimator by analyzing the complete information encoded 
in the time-energy distribution of individual photons in the flare, with the aim of 
probing possible systematic energy-dependent time lags induced by QG vacuum re- 
fraction during photon propagation to the Earth, or intrinsic to the source. The true 
shape of the time profile at the source is not known, so we choose the following 
analysis strategy. In general, the short pulse structure of any flare would be blurred 
by an energy-dependent effect on photon propagation. Conversely, one may correct 
for the effects of any given parametric model of photon dispersion, e.g., the linear 
or quadratic vacuum refractive index, by applying to each photon the appropriate 
time shift [15] corresponding to its propagation in a spatially-flat expanding uni- 

z 

verse: At(E) = Hq^(E/Mqgi) J(l + z)h~^(z)dz or similarly for the quadratic case, 



where Ho is the Hubble expansion rate and h(z) = V^a + ^m(1 + z)^. If the cor- 
rect energy-dependent QG shift is applied, the short pulse structure of the emission 
profile is restored. 
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Fig. 1. The ECF from one realization of the MAGIC measure- 
ments with photon energies smeared by Monte Carlo, for the 
case of a vacuum refractive index that is linear in the photon 
energy. 

We implement this analysis strategy in two ways. In one analysis, we consider the 
most active part of the flare, that is distinguished clearly from the uniform back- 
ground, and the QG shift is varied so as to maximize the total energy in this part. 
In the other analysis we use the shape of the flare as extracted from untransformed 
low-energy data. 



2. 1 Energy cost function 



It is well known [26] that a pulse of electromagnetic radiation propagating through 
a linearly-dispersive medium, as postulated here, becomes diluted so that its power 
(the energy per unit time) decreases 0] Any transformation of a signal to reproduce 
the undispersed signal tends to recover the original power of the pulse. If the pa- 
rameter Mqgi, is chosen correctly, the power of the recovered pulse is maximized. 



^ The applicability of classical electrodynamics for estimating the low-energy behavior 
induced by space-time foam [1,2,4-9, 11] and the corresponding pulse-broadening effect 
have been discussed elsewhere (see [12] for details and an explicit example). The dilution 
effects for the linear or quadratic cases may easily be obtained as described in [26, Sect 7.9] 
by applying the dispersion laws co{k) = k[l - ^c/(2Mqgi)] or aj{k) = k[l - k^/iSM^^^)^. 
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Fig. 2. The t/ distribution from fits to tiie ECFs of 1000 realiza- 
tions of the July 9 flare with photon energies smeared by Monte 
Carlo. 

We implement this observation as follows. First, we choose a time interval (tutj) 
containing the most active part of the flare, as determined using a Kolmogorov- 
Smimov (KS) statistic [27]. The KS statistic is calculated from the difference be- 
tween the cumulative distribution function estimated from the unbinned data and 
that of a uniform distribution. The interval ^2) covers the time range where the 
value of the KS difference varies from its maximum over the whole time support 
of the signal down to a negligible value. This procedure determines the proper 
time-width ?2 - h of the most active (transient) part of the flare[i| Having chosen 
this window, we scan over the whole support the time-distribution of all photons 
shifted by At{E) and sum up the energies of photons in the window. For conve- 
nience, we re-parametrize the time shift as At = +tiE or At = +TqE^ respectively, 
with T/ and Tq in s/GeV and s/GeV^ units. The transformation is repeated for many 
values of r/ and Tq, chosen so that the shifts At match the precision of the arrival- 
time measurements, and for each t/ or Tq the scan is performed and the maximal 
summed energy in a window of width ?2 - ?i is obtained. The maximal energies as 
a function of t/ or Tq define the 'energy cost function' (ECF). The position of the 
maximum of the ECF indicates the value of t/ or Tq that best recovers the signal, 
in the sense of maximizing its powerQ This procedure is applied to 1000 Monte 

^ The time interval chosen agrees very well with the spread of a Gaussian fit to the profile 
of the binned data, as well as the more complicated profile used in [23]. 
^ Varying slightly the boundaries of the interval {t\ ; t2) has a negligible effect on the posi- 
tion of the maximum. We take into account the difference between the width at the Earth 
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Carlo (MC) data samples generated by applying to the measured photon energies 
the (energy-dependent) Gaussian measurement errors. 

Fig. \T\ shows the ECF for one such energy-smeared MC sample. It exhibits a clear 
maximum, whose position may be estimated by fitting it with a Gaussian profile 
in the peak vicinity. Fig. [2] shows the results of such fits to the ECFs with t/ for 
the 1000 energy- smeared realizations of the July 9 flare. From this distribution we 
derive the valuer/ = (0.030+0.012) s/GeV, where Mqgi = 1.445x10^'' s/t/, leading 
to a lower limit Mqgi > 0.21 x 10^^ GeV at the 95% C.l|Z] The same procedure 
applied to the ECF obtained using leads to = (3.71 + 2.57) x 10"'' s/GeV^ 
where Mqgi = 1-222 x 10^ (s/T^)'^^ corresponding to Mqgi > 0.26 x 10" GeV at 
the 95% C.L. While our results for the June 30 flare have similar sensitivities and 
are compatible, they cannot be used to strengthen our results, as this flare is not 
very significant. 

2.2 Likelihood function 

We have confirmed this result using another technique to study the energy-dependent 
delay signal in the data. It is motivated by the initial time and energy-binned analy- 
sis performed in [23], which we used to check that the light-curve is well described 
by a simple Gaussian profile, superimposed on a time-independent background. We 
compute a likelihood function X based on the probability of a photon to be observed 
with energy E and arrival time t, using variables describing the energy spectrum at 
the source, the time distribution at emission obtained from the measured arrival 
times of the photons assuming an adjustable energy-dependent propagation delay, 
and the energy resolution of the detector, which is modelled as a Gaussian [28]. To 
describe the photon energy at the source a simple power law T{E^) ~ is taken, 
with yS = 2.7 for the time-uniform part of the flare and 2.4 for the flaring part [23]. 

The likelihood function is fitted to the July 9 data minimizing - log X as a func- 
tion of four parameters: (i) the energy-dependent delay parameterized in terms of 
Mp/Mqgi, (ii) the position of the intrinsic maximum of the Gaussian flare, (iii) 
its width and (iv) the normalization of the time-independent background compo- 
nent in arbitrary units. The best fit yields M^/Mqgi = 8.2;^34, corresponding to 
T/ = (0.048 + 0.021) s/GeV. The shape of the function = -2 log X + const around 
the minimum in these variables is quite parabolic almost up to the 2-cr level. In view 
of the correlations with these parameters, the sensitivity to ti would be improved if 
they were known more precisely. 

A similar procedure in the case of a quadratic dependency gives = (4.60+5.46)x 
and at the source, also negligible. 

^ We propagate the large errors by using the ±l-o" range of Mp/Mqgu to estimate the 
+l-cr range of Mqcu- 
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Fig. 3. Comparison of the MAGIC measured lightcurve at low and high energies 
with the prediction given by the best set of parameters found using the likelihood 
method, and binning the data and the likelihood function in the same manner. 



10-6 s/GeV2. 

Fig. [3] shows that the X function gives a good overall fit to the data: binning in time 
and energy both the data and the X function, we find ;(f^/NDF ~ 1 . 



2.3 Crosscheck with Monte Carlo data 



To check the robustness of the ECF and likelihood analyses, we simulated several 
MC test samples with two components: (a) a time-independent background with the 
same energy spectrum as the measured data before the flare, and (b) a superposed 
signal generated at the source with an energy spectrum similar to that observed 
during the flare and an energy-independent Gaussian time distribution, each with 
the same numbers of photons as in the measurement. We then calculated the ar- 
rival times of all photons using various dispersion models and parameters, taking 
into account the MAGIC energy resolution. For each dispersion model and param- 
eter, we generated 1000 incarnations, using different random seeds. These samples 
were then analyzed blindly, and the encoded effects were recovered successfully 
by the two estimators within the expected uncertainties. In addition, the analysis 
techniques were applied to MC samples with no energy-dependent dispersive sig- 
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nal encoded, and found no effect, and both techniques also returned null results 
when applied to Mkn 501 data from outside a flare. These tests confirm the numer- 
ical sensitivities of the analyses and the estimates of the uncertainties given above. 
For the likelihood method, additional checks have been performed [28] assuming 
different flare energy spectra and shapes, besides the Gaussian one discussed here, 
which also fit reasonably well the binned data (c.f. Fig.©. 



3 Conclusions 

The probability of the zero-delay assumption relative to the one obtained with the 
ECF estimator is P = 0.026. The observed energy-dependent delay thus is a likely 
observation, but does not constitute a statistically firm discovery. The results of 
the two independent analyses of the July 9 flare of Mkn 501 are quite consistent 
within the errors. Their results exhibit a delay between y-rays of different energies, 
Ti = (0.030 + 0.012) s/GeV, corresponding to a lower limit Mqgi > 0.21 x 10'** GeV 
at the 95% C.L. We also find a quadratic delayr^ = (3.71+2.57)xl0"^ s/GeV^, and 
Mqg2 > 0.26 X 10'^ GeV at the 95% C.L., far beyond previous limits on a quadratic 
effect in photon propagation [1 1, 14, 18]. These numbers could turn into a real mea- 
surement of Mqgi,25 if the emission mechanism at the source were understood and 
the observed delays were mainly due to propagation. We cannot exclude, however, 
the possibility that the delay we find, which is significant beyond the 95% C.L., is 
due to some energy-dependent effect at the source0 However, we can exclude the 
possibility that the observed time delay may be due to a conventional QED plasma 
refraction effect induced as photons propagate through the source. This would in- 
duce [29] At = D{a^T^ /6q^) lv?{qT /ml), where a is the fine-structure constant, q is 
the photon momentum, T is the plasma temperature, is the mass of electron, D 
is the size of the plasma, and we use natural units: c,h = 1. Plausible numbers such 
as r ~ 10"^ MeV and D ~ 10^ km (for a review see [30]) yield a negligible effect 
for (7 ~ 1 TeV. Exclusion of other source effects, such as time evolution in the mean 
emitted photon energy, might be possible with the observation of more flares, e.g., 
of different AGNs at varying redshifts. Observations of a single flare cannot distin- 
guish the quantum-gravity scenarios considered here from modified synchrotron- 
self-Compton mechanisms [23,31]. However, this pioneering study demonstrates 
clearly the potential scientific value of an analysis of multiple flares from different 
sources. The most promising candidate for applying the analyses proposed here is 
the flare from PKS 2155-304 detected recently by H.E.S.S. [32]. Unfortunately the 
occurrence of fast flares in AGNs is currently unpredictable, and since no correla- 
tion has yet been established with observations in other energy bands that could be 
used as a trigger signal, only serendipitous detections are currently possible. 



^ Note that if the observed energy-dependent time shift is explained by some source effect, 
the lower limit on Mqq would rise. 
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